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Abstract 

Folding of the triangular lattice in a discrete three-dimensional space is studied numerically. 
Such "discrete folding" was introduced by Bowick and co-workers as a simplified version of the 
polymerized membrane in thermal equilibrium. According to their cluster- variation method (CVM) 
analysis, there appear various types of phases as the bending rigidity K changes in the range 
— oo < K < oo. In this paper, we investigate the K < regime, for which the CVM analysis 
with the single-hexagon-cluster approximation predicts two types of (crumpling) transitions of 
both continuous and discontinuous characters. We diagonalized the transfer matrix for the strip 
widths up to L = 26 with the aid of the density- matrix renormalization group. Thereby, we 
found that discontinuous transitions occur successively at K = —0.76(1) and —0.32(1). Actually, 
these transitions are accompanied with distinct hysteresis effects. On the contrary, the latent-heat 
releases are suppressed considerably as Q = 0.03(2) and 0.04(2) for respective transitions. These 
results indicate that the singularity of crumpling transition can turn into a wea/i-first-order type 
by appreciating the fluctuations beyond a meanfield level. 



I. INTRODUCTION 



A model of "discrete foldinEr" was introduced as a possible lattice realizationfdiscretized 
^version) for the tethered (polymerized) membranes in thermal equilibrium 



The constituent molecules of the tethered membranes are connected with each other 



via polymerization j9 



10 



111. Hence, the in-plane strain is subiected to finite shear moduli. 

nn 

In this sense, the tethered membranes are stiff, compared with the fluid membranes [Ij, 113|| 
for which the constituent molecules are diffusive, and cannot support a shear. Reflecting 
this peculiarity, the tethered membranes get flattened macroscopically at sufficiently low 
temperatures (or equivalently large bending rigidities) Qj. This phase transition is called 
the crumpling transition. (Because the crumpling transition occurs even in the absence of 
the self-avoidance, we neglect this effect throughout this paper.) The flat phase is charac- 
terized by the long-range orientational order of the surface normals. It would be remarkable 
that such an 0(3) symmetry, namely, an orientational order, breaks spontaneously for a 
two-dimensional manifold. Actually, it has been known that the fluid membranes are always 
crumpled irrespective of the temperatures. To clarify the nature of the crumpling transition 
of the tethered membranes, a good deal of theoretical analyses have been reported so far. 
However, there still remain controversies whether the singularity belongs to a continuous 
phase transition [l^ 
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companied with appreciable latent-heat release 
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27| or a discontinuous one ac- 



29|. Actually, in numerical simulations, 
it is not quite obvious to rule out a possibility of weaA;-first-order phase transition 0,0. 
As would be anticipated, an approach via the discrete-folding model 

mm 

IS promis- 
ing to resolve this longstanding issue: In fact, the discrete-folding model reduces to a mere 
Ising magnet ^ with a dual transformation; we explain the details afterward. The Ising 
variables, however, are subjected to a constraint, which makes the thermodynamics rather 
nontrivial. Nevertheless, once formulating the discrete folding in terms of the Ising variable, 
we are able to resort to a variety of techniques developed in the study of magnetism: Cirillo 
and co-workers carried out an extensive cluster-variation method (CVM) analysis with the 

riQ 

single-hexagon-cluster approximation p, Remarkably enough, according to their CVM 
solution, various types of phases emerge as the bending rigidity K changes; see Fig. ?? for 
a schematic drawing of the phase diagram. As depicted in the figure, the symbols, Kp^, 
Kg3, and Kc3, specify the respective crumpling-transition points separating the piled-up. 
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tetrahedral, octahedral, and totally- flat phases; we followed the notation in Ref. jsf- 

Stimulated by their result, in an earlier paper we performed a first-principles sim- 
ulation for the discrete folding. Our simulation scheme is based on the transfer-matrix 
diagonalization with the aid of the density-matrix renormalization group 33ll34l.l35ll36l.l37|. 
In Ref. |32], we dwelt on the > regime, and investigated the crumpling transition at 
K = Kc3'- We estimated the transition point as Kc3 = 0.195(2), and estimated the amount 
of the latent heat as Q = 0.365(5). Moreover, we confirmed that for K > Kc3, the en- 
tropy disappears completely as predicted by the CVM analysis ^. That is, the membrane 
stretches out into a totally fiat sheet in that regime. In this sense, the transition at K = Kcs 
is rather peculiar, and it would be less relevant to that of the realistic tethered membranes. 
Beside slight quantitative discrepancies concerning the folding entropy at = and the 
amount of the latent heat, our simulation data support the CVM predictions, implying that 
the membrane undulations beyond a single-hexagon cluster are not quite significant in the 
K > regime. 

In this paper, we devote ourselves to the K < regime, where the CVM analysis predicts 
two types of crumpling transitions [8j. Our aim is to examine the singularities of these 
transitions with the density-matrix renormalization group for L < 26: It is conceivable that 
the order of the transition changes or that even the transition disappears by appreciating the 
fluctuations beyond the single-hexagon-cluster level. As a consequence, we found that the 
crumpling transitions occur at K = —0.76(1) and —0.32(1), and the singularities are both 
discontinuous. Actually, we estimated the amounts of the latent heat as Q = 0.03(2) and 
0.04(2), respectively. Our result indicates that the singularity of the crumpling transition 
can change into a (weak) first-order type by fully appreciating the fluctuations. 

In fairness, it has to be mentioned that the discrete folding has been studied extensively 
with computer simulations other than the density-matrix renormalization group, namely, 
the conventional full-diagonalization scheme 0, 1^ and the Monte Carlo method JssJ: In 
the full-diagonalization calculation, the tractable systems sizes are limited within L < 6, 
and are not feasible for identifying the characters of the transitions. In fact, as for the 
transition at K = Kp^, it was speculated that no sign of the transition could be captured 
until L = 8, which is far beyond the ability of the full-diagonalization scheme. On the 
other hand, the Monte Carlo simulation suffers from a very slow relaxation to the thermal 
equilibrium (glassy behavior). Although such a metastability is a fascinating topic in its 
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own right, it brings about a severe obstacle to an efficient sampling in thermal equilibrium. 
(Such slow relaxation is to be attributed to the constraint to the Ising variables.) We stress 
that the density-matrix renormalization group is capable of treating large system sizes, and 
free from the slow-relaxation problem. 

The rest of this paper is organized as follows. In the following section, we give an account 
of the discrete folding, following the presentation in Ref. Q]. Then, we explain the density- 
matrix renormalization group adopted to the discrete-folding model 3^. In Sec. ??, we 
present the numerical results. The last section is devoted to the summary and discussions. 



II. DIAGONALIZATION OF THE TRANSFER MATRIX WITH THE DENSITY- 
MATRIX RENORMALIZATION GROUP: A BRIEF REMINDER OF REF. 



In this section, we explain an outline of our simulation scheme. Full account of details 
will be found in an earlier paper 0|. To begin with, we introduce the discrete folding of 
the triangular lattice in a discretized three-dimensional space Q]. 

n 

A. Discrete folding in three dimensions ||^] 

Let us explain how the triangular lattice is folded (embedded) in a discrete three- 
dimensional space, namely, the face-centered cubic (fee) lattice. As shown in Fig. 3 of 
Ref. P], the fee lattice is viewed as a close packing of two types of polygons, namely, the 
octahedrons and tetrahedrons, whose vertexes are shearing the fcc-lattice points; the faces 
of these polygons are all equilateral triangles. Hence, wrapping these polygons with a sheet 
of the triangular lattice, we are able to embed the sheet arbitrarily in the fee lattice. More 
specifically, the relative fold angles 9 between adjacent triangles are taken from four possi- 
bilities, namely, "no fold" {6 = tt), "complete fold" {9 = 0) "acute fold" {9 = acos{l/3)), 
and "obtuse fold" {9 = acos{-l/3)). 

The above discretization leads to an Ising-spin representation for the discrete folding. 
An efficient representation, the so-called gauge rule, reads as follows 5]: We place two 
types of Ising variables, namely, cxj and Zi, at each triangle (rather than each joint); see 
Fig. ?? (a). The gauge rule states that these Ising spins specify the joint angle between 
the adjacent triangles. That is, provided that the z spins are antiparallel (,21^2 = —1) for 
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a pair of adjacent neighbors, the joint angle is either acute or obtuse fold. Similarly, if 
o"iO"2 = —1 holds, the relative angle is either complete or obtuse fold. Note that the above 
rules specify the joint angle unambiguously. To summarize, we introduced Ising variables 
placed at each triangle. Because the structure dual to the triangular lattice is hexagonal, the 
Ising variables constitute an Ising model on the hexagon lattice. Hence, the transfer-matrix 
strip looks like that drawn in Fig. ?? (b). The row-to-row statistical weight T^ai,zi}.{a',z'} 
yields the transfer-matrix element. However, according to Ref. |5j, the Ising variables are 
not quite independent, but subjected to some constraints for each hexagon. The transfer- 
matrix element is given by the following form with extra factors that enforce the constraints 

S; 

T{z,.U.',.'} = (n U,V,)exp{-H/T), (1) 

with, 

Uj = 6{a2j-2 + o-2j-i + (T2j + (J2j_i + o-2i + ^2j+l 3, 0), (2) 

and, 

2 

Vj = n SMz2j, Z2j-i, Z2j-2, 4i-i> 4j> ^od 2, 0). (3) 

c=l 

Here, S{m,n) denotes Kronecker's symbol, and ac is given by, 

6 i 

ac{zi, . . . , ze) = XI 9 (1 ~ ZiZi+i)S{co + '^ai-c mod 3, 0). (4) 

i=i ^ j=i 

The Boltzmann factor exp{—H/T) in Eq. (??) is due to the bending-energy cost. Hereafter, 
we choose the temperature T as the unit of energy; namely, we set T = 1. As usual, the 
bending energy is given by the inner product, cos 6ij of the surface normals of adjacent 
triangles. Hence, the bending energy is given by the formula; 

H = -0.5 cos 9ij = -0.5 Y,^Kai(T j{l + 2ziZj), (5) 

with the bending rigidity K. Here, the summation Y.{ij) runs over all possible nearest- 
neighbor pairs (ij) around each hexagon. (The overall factor 0.5 is intended to reconcile the 
double counting.) The above completes a prescription to construct the transfer matrix. In 
the following, we explain how we diagonalized the transfer matrix. 
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B. Diagonalization of the transfer matrix for the discrete folding: An application 
of the density-matrix renormalization group js^ 

As would be apparent from Fig. ?? (b), a unit cell contains four Ising variables. Therefore, 
the size of the transfer matrix grows rapidly in the form 16^ x 16^ with the system size L. 
Hence, the diagonalization of the transfer matrix requires huge computer-memory space. 
(Actually, with the conventional full-diagonalization scheme, the tractable systems sizes are 
limited within L < 6 .) In order to cope with this difficulty, in an earlier paper |32| |. 

we developed an alternative (memory saving) diagonalization method based on the density- 



matrix renormalization group. 



33, 



gorithm, 
3911 for a 



34L l35| . In the following, we outline the a 
placing an emphasis on the changes specific to this problem; refer to a textbook 
pedagogical guide. 

The density-matrix renormalization group is based on the idea of the (computer aided) 
real-space decimation. We presented a schematic drawing of one operation of the real- 
space-decimation procedure in Fig. ??. Through the decimation, the block states and the 
adjacent spin variables (hexagon) are renormalized altogether into a new block states. Note 
that through the decimation, the number of states for the (renormalized) block is truncated 
within m; the parameter m sets the simulation precision. Hence, sequential applications of 
the procedure enable us to reach very long system sizes. In order to reduce the truncation 
error, we need to retain m significant bases. According to the criterion advocated in Ref. 
Q, ... ,„ .a., a. own fro. eige.ta.e. ot the (loca„ den... ,.a... w.t. 
dominant statistical weights (large eigenvalues) {wa} {en: integer index) |39l |. 

This is a good position to address a few remarks regarding the changes specific to the 
present problem. First, in our simulation, in order to reduce the truncation error of the 



real-space decimation, we adopted the "finite-size method" |39|. Secondly, we applied 
a magnetic field h at an end of the transfer-matrix strip. Namely, we incorporated an 
additional Hamiltonian of the form —h6cr,iSz,i with h = 0.05 and the edge spins {a, z). This 
trick was utilized in the preceding full-diagonalization calculation j^, and it aims to split 
off the degeneracy caused by the trivial four-fold symmetry (cr = ±l,2; = ±l)of the overall 
membrane orientation (gauge redundancy). 
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III. NUMERICAL RESULTS 



In this section, based on the algorithm explained above, we perform computer simula- 
tions for the discrete-folding model in the negative bending rigidity {K < 0) regime. As a 
preliminary survey, we calculate the order parameters, Eqs. (??)-(??), for a wide range of 
K. Thereby, we examine the phase diagram predicted by the CVM analysis as depicted in 
Fig. ??. We will also provide a distribution of the density-matrix eigenvalues in order to 
elucidate the reliability of the present scheme. These preliminary analyses are followed by 
large-scale simulations to determine the singularities of the crumpling transitions precisely. 



A. Preliminary survey: A check of the rehabihty of the simulation scheme 

In Fig. ??, we plotted the following order parameters (local magnetizations) for the 
system size L = 14 and the number of states remained for a renormalized block m = 12; 

O = (a,) (6) 
T = (zas) (7) 
P = {z), (8) 

where the bracket (■ ■ ■) denotes the thermal average, and the symbol as stands for the stag- 
gered component of the moment a. (Non-zero magnetizations are induced by the symmetry- 
breaking field at an end; see Sec. |Tllfor details.) As depicted in Fig. ??, the CVM analysis 
predicts successive crumpling transitions with varying the bending rigidity K, and these 
phases are characterized by the above order parameters. Indeed, our simulation data support 
this claim; namely, there appear crumpling transitions at i^' ^ —0.9 and —0.3 successively. 
However, on closer inspection, our simulation data exhibit clear evidences of the hysteresis 
effects for both transitions. Hence, contrary to the CVM scenario, both transitions should 
be of discontinuous characters. 

Before going into detailed analysis of the hysteresis effects, let us check the reliability of 
the scheme. In Fig. ??, we show the distribution of the statistical weights (eigenvalues of 
the density matrix) {wa} for K = —0.2, L = 15 and m = 15. We see that the statistical 
weight drops very rapidly. In other words, a majority of the statistical weight concentrates 
on a few dominant (relevant) states. 
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From the data, we are able to read off a reliability of the numerical simulation: Because 
we discard those states with the indexes a > m through an operation of the renormalization- 
group procedure, the statistical weight Wq, at a = m yields an indicator of the truncation 
error. To be specific, the fifteenth state, for instance, exhibits a very tiny statistical weight 
5 X 10~^. This truncation error accounts for a precision of the free energy. The 
reliability of other quantities such as the internal energy E could be enhanced by a factor 
of 10; namely, the precision would be of the order ~ 5 x 10~^. The data scatter of Figs. 
?? and ??, for instance, should be attributed to this uncertainty. Moreover, it is to be 
noted that the strength of bending rigidity set in Fig. ??, namely, \K\ = 0.2, is rather 
small. Therefore, upon stiffening the rigidity li^l further, the fluctuations get suppressed, 
and correspondingly, the truncation error improves. In particular, for \K\ > 0.7(~ IfCpsI), 
such an improvement appears to be pronounced. Encouraged by these findings, we proceed 
to large-scale simulations to determine the singularities of the crumpling transitions. 

B. Singularity of the crumpling transition at K = Kg3 

In the above, we observed clear evidences of the onsets of the crumpling transitions at 
K ^ —0.9 and —0.3. Contrary to the CVM scenario, these transitions are both discontin- 
uous. Here, performing large-scale simulations, we confirm that the transitions are indeed 
discontinuous ones accompanied with appreciable latent-heat releases. 

We first consider the crumpling transition at K = Kgs, namely, the transition separating 
the tetrahedral and octahedral phases. In Fig. ??, we plotted the free energy per triangle 
/_^^^ for L = 20, and m = 15. (See our preceding paper [32| for an algorithm to calculate the 
free energy reliably.) As indicated in the plot, we swept the parameter K in both descending 
and ascending directions. Correspondingly, the indexes of and specify the parameter- 
sweep directions. From the data, we observe a clear hysteresis effect; the intersection point 
of and /<_ yields the location of the first-order phase transition. Therefore, we confirm 
that the crumpling transition at = Kg3 is indeed a first-order phase transition. 

In order to determine the transition point more precisely, in Fig. ??, we plotted the 
excess free energy — for various values of L and m. From the figure, we see that the 
data for L > 20 converge to the thermodynamic limit satisfactorily. From the coexistence 
condition /<_ — = 0, we estimate the transition point as = —0.32(1). The location of 
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the transition point is in good agreement with the CVM result K = —0.294. 

In order to estimate the latent heat, in Fig. ??, we plotted the residual internal energy 
— per triangle for various values of L and m. Similar to the above, the data, E^ 
and E^, are calculated by sweeping K in ascending and descending directions, respectively. 
The residual internal energy at the transition point K = —0.32(1) yields the latent heat. 
The data for L > 20 appear to reach the thermodynamic limit satisfactorily. From these 
data, we determine the latent heat a.s Q = 0.04(2). The present estimate is not consistent 
with the CVM result Q = 0.14; this value is taken from Fig. 12 of Ref. ^. We see that 
the discrepancy between the present estimate and the CVM analysis with a single-hexagon- 
cluster approximation is rather severe, and it may indicate that the fluctuations beyond a 
hexagon cluster are not quite negligible. Such fluctuations contribute to smearing out the 
latent-heat release significantly. 

Although the latent-heat release acquires such a notable suppression, there appears a 
distinct hysteresis effect as shown in Figs. ?? and ??. This_^eculiarity may be reflected 
in the preceding fuU-diagonalization calculation for L < 6 [8]: This calculation captures 
an evidence of the onset of the transition by a steep increase of the order parameter T, 
whereas the character of the singularity was remained unsolved. Actually, such a small 
amount of latent heat could hardly be resolved by conventional approaches. Nevertheless, 
such peculiarities, namely, a small latent heat and a distinct hysteresis effect, suggest that 
the barrier between the coexisting states is hardly surmountable, albeit these states are 
rather close to each other in the configuration space. To the best of our knowledge, such 
a feature is specific to the crumpling transition, and it may be related to the anomalous 
metastability (glassy behavior) observed in the Monte Carlo simulation [ssj ]. 

Lastly, let us address a technical remark on the hysteresis effect. The hysteresis cannot 
be reproduced by the conventional fuU-diagonalization scheme. In the density-matrix renor- 
malization group, on the contrary, in the course of simulation, past information is encoded 
(retained) in the renormalized block through the sequential applications of normalizations. 
Because of this "memory effect" 0|, the density-matrix renormalization group reproduces 
a hysteresis behavior. Beside this advantage, the method admits reliable estimate for the 
free energy. These advantages are significant to determine the first-order-phase-transition 
point reliably. 
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C. Singularity of the crumpling transition at K = Kp^ 



Let us turn to surveying the singularity at K = Kps, namely, the crumpling transition 
separating the piled-up and tetragonal phases; note that according to CVM, the transition 
should be continuous. 

In Fig. ??, we plotted the excess free energy — for various values of L and m. 
The data exhibit a clear hysteresis effect. Hence, we confirm that contrary to the CVM 
scenario, the transition is of a discontinuous character. From the coexistence condition 
— = 0, we determine the transition point as i^' = —0.76(1). The location of the 
transition point is consistent with the CVM estimate K = —0.852. Hence, apart from the 
singularity characterization, the CVM treatment appears to be justified. 

The discrepancy between the present numerical result and the CVM analysis suggests 
that just like the above mentioned case of = Kqs, the fluctuations beyond a hexagon 
cluster are important. As a matter of fact, the preceding fuU-diagonalization analysis for 
L < 6 could not even detect a sign of singularity around K = Kp^, and the authors 
speculated that the singularity would not be captured correctly until L = 8. In other words, 
we need to consider long-length-scale fluctuations exceeding L = 8, at least, in order to 
capture the characteristics of this transition. In this sense, the crumpling transitions in the 
K < regime are much more subtle than that of > studied previously in Refs. 3, ■ 

Moreover, as shown in Fig. ??, the behavior of the excess free energy — is asym- 
metric with respect to the transition point K = —0.76(1); namely, the absolute value of the 
excess free energy is enhanced in the K < —0.76 side, whereas in the other size K > —0.76, 
it remains suppressed. (Note that as for K ~ Kg3, for both sides of the transition, the excess 
free energy behaves quite similarly (symmetrically) as shown in Fig. ??.) A large amount 
of the excess free energy reflects an existence of strong metastability. In fact, as shown in 
Fig. ??, the order parameters exhibit notable hysteresis effects in the regime K < —0.76. 
That is the reason why a symptom of the transition implied by the hysteresis curve deviates 
significantly from the true transition point determined by the coexistence condition (zero 
excess free energy) f^ — f^ = 0. We suspect that the hexagon-cluster- approximation CVM 
estimate Kp^ = —0.852, which is slightly smaller than ours Kp^ = —0.76(1), could also be 
biased by this asymmetry. Nevertheless, as mentioned in the Monte Carlo simulation js^, 
such a metastability, reminiscent of that of the spin glasses, is quite peculiar, and further 
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details have be clarified so as to confirm the present simulation result. 

In order to estimate the latent heat, we plotted the residual internal energy — 
in Fig. ??. From the residual internal energy at the transition point, we estimate the 
latent heat a.s Q = 0.03(2). As noted above, such a finite latent-heat release is due to the 
fluctuation effect beyond a meanfield level. It is to be noted however, in the present case 
of K = Kps, the fluctuation effect leads to enhancing the latent-heat release rather than 
smearing it out as in the above mentioned case of i^' = Kgs- 

We notice that the amount of the latent heat is considerably suppressed. (Note that the 
latent heat at K = Kc3{> 0) is much larger Q = 0.365(5) 32j. ) Our simulation result 
suggests that the long-length-scale fluctuations beyond the meanfield treatment can drive 
the transition to a ti;eaA;-first-order type. As mentioned in Introduction, the singularity of 
the crumpling transition for (realistic) membranes is arousing controversies. In this respect, 
the present result cautions us not to rule out the possibility of weak-first-order type very 
easily so as to avoid misidentifying it with a continuous one of an unknown universality class. 
Nevertheless, we suspect that the transition in the i^' < side is relevant to the crumpling 
transition of (realistic) tethered membranes rather than the transition in the other side 
K < 0, whose discontinuous character is too pronounced. 

n 

According to the analytical theory [8|] , the crumpling transition at K = Kp^ is equivalent 
to that of the two-dimensional (planar) folding H S which exhibits a continuous 
transition in < JJ]. Hence, provided that this equivalence holds, the crumpling transition 
at = Kp3 should be continuous as well. This mapping is based on the idea that for large 
li^'l, the complete and acute folds are dominant, and within this restricted configuration 
space, the discrete folding is equivalent to the planar folding. However, our simulation 
result does not support this claim, and rather, it suggests that the restricted-configuration 
picture would not be validated. Possibly, a slight violation of such a restriction alters the 
transition singularity into a wea/c-first-order type. In other words, in reality, various types 
of fluctuations are promoted around K = Kp^, and such fluctuations give rise to a finite 
latent-heat release. 
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IV. SUMMARY AND DISCUSSIONS 



We investigated the singularities of the crumphng transitions of the discrete-folding model 
in the negative-bending-rigidity regime i^' < 0. We adopted the density-matrix renormal- 
ization group to diagonalize the transfer matrix for the strip widths up to L = 26; the 
system sizes treated here are substantially larger than those tractable with the conventional 
fuU-diagonalization scheme L < 6 [8]. Actually, in Ref. ^, it was speculated that the 
singularity at K = Kp^ would not be captured correctly until L = 8. Hence, it is significant 
to consider long- length-scale undulations exceeding L = 8. In order to check the reliability 
of the present simulation, we calculated the distribution of the density-matrix eigenvalues 
(statistical weights) for L = 15 and K = —0.2. From the data, we found that the truncation 
error due to the renormalization group (real-space decimation) procedure is of the order 
~ 5 X 10~^ for m = 15 (number of bases remained for a renormalized block). This result is 
encouraging in the sense that the truncation error is smaller than the finite-size corrections, 
and in practice, almost negligible. 

As a preliminary survey, we calculated the order parameters, Eqs. (??)-(??), for a wide 
range of K{< 0). The overall behaviors support the CVM scenario as depicted in Fig. ??; 
namely, there appear two crumpling transitions around K ^ —0.9 and —0.3, which separate 
the piled-up, tetragonal and octahedral phases. However, on closer inspection, we observe 
clear hysteresis effects for both transitions. That is, contrary to the CVM scenario, both 
singularities should be discontinuous. 

Performing large-scale simulations, we determined the transition points as = —0.76(1) 
and —0.32(1). At respective transition points, we estimated the amounts of the latent 
heat as Q = 0.03(2) and 0.04(2). Comparing the present results with the CVM estimates 
with the single-hexagon-cluster approximation, namely, Q = {Kp^, = —0.852) and 0.14 
{Kg3 = —0.294), we notice that the discrepancies between them are rather conspicuous. 
The present simulation results suggest that the fluctuations beyond a hexagon cluster are 
significant. According to Ref. [q], the transition at K = Kp^ should be a continuous 
one, provided that complete and acute folds are dominant around the transition point; 
within this restricted configuration space, the discrete folding is equivalent to the planar 
(two-dimensional) folding, which exhibits a continuous transition in i^' < 0. However, 
our simulation does not support this claim, because the transition turns out to be of a 
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discontinuous character. 

In an earlier paper js^, we studied the crumphng transition in K > 0, and reported that 
it belongs to a first-order phase transition with the latent heat Q = 0.365(5). Consequently, 
it turns out that all crumpling transitions of the discrete-folding model are discontinuous. 
However, it is notable that the amounts of the latent heat in K < are much smaller than 
that of the other side K > 0. In other words, the singularities of the crumpling transitions 
in i^' < are all wea/c-first-order phase transitions. As mentioned in Introduction, on 
the contrary, a majority of Monte Carlo simulation for (realistic) tethered membranes has 
reported that the crump ling tra nsition should be a second-order phase transition [3, 1^ 



17 



21 



22 



23,^ 
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26 



27l |. However, as noted in Refs. j30|,|3l|, it is not quite 



obvious to rule out the possibility of weak first-order transition in practice. Hence, it would 
be desirable to set forward the present analysis to other discretized-membrane models, such 
as the square-diagonal- lattice folding and its variants 4^ 3,13' order to see whether the 
crumpling transitions are still of discontinuous characters (weak first order); the discretized- 
membrane model admits detailed analysis on the character of the crumpling transition. This 
problem will be addressed in future study. 
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FIG. 1: (a) A p: 
method (CVM) 



)hase diagram of the discrete folding model determined with the cluster- variation 
~~\. Various phases appear as the bending rigidity K changes, (b) A schematic 
drawing of the behaviors of the order parameters, Eqs. (??)-(??), determined with CVM _8]; this 
result is to be compared with our first-principles-simulation data as shown in Fig. ??. Note that 
according to CVM, the transition at Kp^ (Kgs) should be a continuous (discontinuous) one. 
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FIG. 2: (a) We consider a discrete folding of the triangular lattice embedded in the three- 
dimensional space. In order to specify the fold angle, we place two types of Ising variables such as 
Zi and (Tj at each triangle rather than at each joint (gauge rule ^). Hence, hereafter, we consider 
a spin model defined on the dual (hexagonal) lattice, (b) A construction of the transfer matrix. 
The row-to-row statistical weight yields the transfer-matrix element. The explicit formula is given 
by Eq. (??). The transfer matrix is diagonalized with the density-matrix renormalization group; 
see Fig. ??. 
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Blach states 1--m 




FIG. 3: Schematic drawing of the density-matrix renormahzation group (DMRG) procedure. 
Prom the drawing, we see that through one operation of DMRG, a "block" and the adjacent sites 
(hexagon) arc rcnormalizcd altogether into a new renormalized "block." At this time, the number 
of block states is retained within m; see text for details. In this manner, we can diagonalize a 
large-scale transfer matrix through sequential applications of DMRG. 
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FIG. 4: The order parameters, Eqs. (??)-(??), are calculated with the density-matrix renormal- 
ization group for the system size L = 14 and the number of states kept for a block m = 12. The 
arrows indicate the directions of the parameter sweeps. We notice that the behaviors are similar 
to those obtained with the CVM analysis (see Fig. ??). However, contrary to the prediction, the 
simulation data exhibit clear evidences of the hysteresis effects for both transitions, implying that 
these transitions are both discontinuous. 
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FIG. 5: Distribution of the eigenvalues (statistical weights) {wa} (a, integer index) of the (local) 
density matrix Q, y, Q is shown. The simulation parameters are chosen from the bending 
rigidity K = —0.2, the number of states kept for a block m = 15, and the system size L = 15. We 
see that Wa drops very rapidly for large a. Because we discard the irrelevant states a > m, the 
weight Wa at a = m yields an indicator of the truncation error. We attain a precision of the order 
5 X 10^"^ under setting m = 15 for instance. 
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FIG. 6: Free energy per triangle f—>,^ is plotted for L = 20 and m = 15. These subscripts, 
and specify the directions of the parameter sweeps, namely, the ascending and descending 
directions, respectively. We see that at ~ —0.3, a first-order phase transition occurs. 
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FIG. 7: Excess free energy /«_ - is plotted for (+) L = 18 and m = 15, (x) L = 20 and 
m = 15, (*) L = 20 and m = 20, and (□) L = 22 and m = 15. Prom the coexistence criterion 
f^ — f^ = 0, we determined the transition point as K = —0.32(1). 
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FIG. 8: Residual internal energy — per triangle is plotted for (+) L = 18 and m = 15, 
(x) L = 20 and m = 15, (*) L = 22 and m = 15, and (□) L = 24 and m = 15. From the residual 
internal energy at the transition point, we estimate the latent heat as Q = 0.04(2). 
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FIG. 9: Excess free energy /<_ - /_» is plotted for (+) L = 18 and m = 20, (x) L = 20 and 
m = 15, (*) L = 22 and m = 15, and (□) L = 22 and m = 20. Prom the coexistence criterion 
/<_ — = 0, we determined the transition point as K = —0.76(1). 
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FIG. 10: Residual internal energy — per triangle is plotted for (+) L = 19 and m = 15, 
(x) L = 20 and m = 15, (*) L = 22 and m = 20, and (□) L = 26 and m = 15. From the residual 
internal energy at the transition point, we estimate the latent heat as Q = 0.03(2). 
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